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Objectives 

Our objective is to help explain how the earliest ancestors of contemporary cells (protocells) 
performed their essential functions employing only the molecules available in the protobi- 
ological milieu. Our hypothesis is that vesicles, built of amphiphilic, membrane-forming 
materials, emerged early in protobiological evolution and served as precursors to protocells. 
We further assume that the cellular functions associated with contemporary membranes, 
such as capturing and transducing of energy, signaling, or sequestering organic molecules 
and ions, evolved in these membrane environments. An alternative hypothesis is that these 
functions evolved in different environments and were incorporated into membrane-bound 
structures at some later stage of evolution. 

We focus on the application of the fundamental principles of physics and chemistry to 
determine how they apply to the formation of a primitive, functional cell. Rather than at- 
tempting to develop specific models for cellular functions and to identify the origin of the 
molecules which perform these functions, our goal is to define the structural and energetic 
conditions that any successful model must fulfill, therefore providing physico-chemical bound- 
aries for these models. We do this by carrying out large-scale, molecular level computer 
simulations on systems of interest. 

Specific systems that we have investigated are: 

1. Simulations of water-membrane systems to yield an accurate description of the elec- 
trical properties of the membranes. This was investigated by calculating the surface 
potentials of aqueous interfaces and the interaction of a number of small polar and 
nonpolar solute molecules with several membrane interfaces. 

2. Simulations of ion transport across a water-membrane interface to investigate how 
“thinning” defects in the membrane effectively increase the permeability of the mem- 
brane by several orders of magnitude over predictions based on simple, dielectric con- 
tinuum models. 
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3. Simulations of small peptides at the water-hexane interface to investigate how the 
interface affects the structure of the peptide, and in particular if short peptides, which 
are disordered in water, can fold into ordered structures at the interface. 


Summary 

Structure of membranes and their permeability to small molecules and ions. Our 

earlier studies on membranes composed of glycerol 1-monooleate (GMO), which are suffi- 
ciently simple to be good models for protobiological membranes, were extended to phospho- 
lipid bilayers composed of dipalmitoylphosphatidylcholine (DPPC) and 1-palmitoyl 2-oleoyl 
sn-glycero 3-phosphocholine (POPC) All these studies show that water readily penetrates 
the head groups but not the hydrophobic, highly disordered core of the bilayer. The surface 
of the bilayer undergoes considerable fluctuations which lead to the formation of thinning 
defects in the membrane. Similar defects provide pathways for unassisted ion transport 
across membranes. Simple monovalent cations, such as Na + , move across the membrane 
well-solvated by water and lipid head groups. The formation of defects increases membrane 
permeability to ions by 13-15 orders of magnitude, rendering the transport feasible even in 
the absence of ion channels or carriers. 

The water-membrane interface was shown to have the unique property that it concen- 
trates many small solutes that exhibit some degree of polarity but are not necessarily am- 
phiphilic. A general theory of this phenomenon was developed based on the balance between 
electrostatic and non-electrostatic contributions to the free energy of transferring these so- 
lutes across the interface. Explaining how small molecules are distributed in membranes has 
broad implications for understanding their permeation across protocellular walls, primitive, 
heterogeneous catalysis and simple, photo-induced charge transfer. 

Peptides at membrane interfaces. In analogy with modern proteins, which fold into 
well-defined structures, we expect that peptides must adopt ordered structures before they 
can perform various protocellular functions. The ability of small peptides to organize at 
aqueous interfaces was examined by investigating peptides composed of two amino acids, 
nonpolar leucine and polar glutamine, but having different sizes and sequences of amino 
acids. Among them were dipeptides, heptamers designed to maximize the interfacial stability 
of an a-helix and a /?-strand, and an undecamer composed entirely of leucine residues. Using 
the undecamer of leucine as a model, the complete folding of a peptide from a disordered 
structure to an a-helix was accomplished for the first time in computer simulations that 
correctly accounted for solvent effects. These simulations revealed principles governing the 
organization of peptides at interfaces: 

1. Short peptides tend to accumulate at interfaces and acquire ordered structures, pro- 
viding that they have a proper sequence of polar and nonpolar amino acids. The 
specific identity of amino acids is less important, a desirable protobiological property. 
Once folded, the peptides may organize into structures suitable for polymerization or 
catalytic activity. The interfacial structure of peptides is mainly determined by the 
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hydrophobic effect, which is manifested at aqueous interfaces as a tendency for polar 
and nonpolar groups of the solute to segregate into the aqueous and nonpolar phases, 
respectively. 

2. If peptides consist of nonpolar residues only, they partition into the nonpolar phase 
and simultaneously fold into an a-helix. Once in the nonpolar environment, they 
can adjust their structure and orientation to changing external conditions. This may 
have provided a simple mechanism of transmitting signals from the environment to the 
interior of a protocell. 


Introduction 

An important step in the origin of life was the organization of material in the prebiotic 
environment into a protocell, which is a hypothetical precursor of modern cells. A protocell 
would have been able to sequester material from the external environment, to carry out 
the chemical reactions necessary to maintain its integrity and allow growth, and to divide 
(thereby reproducing), possibly in the absence of a genome. Since all known cells have 
membrane-like barriers separating the interior from the external environment and since all 
known cells use proteins to catalyze the chemical reactions along their metabolic pathways, 
we envision a protocell that is formed from a simple, lipid bilayer membrane that uses 
peptides short amino acid polymers that are simpler than modern proteins — as primitive 
catalysts. This model of a protocell forms the basis of the research described herein: how 
a bilayer membrane, in which hydrophobic and hydrophilic media exist in close proximity, 
affects the concentration of material in its vicinity and the transport of material across the 
membrane; and how the membrane can affect the structure of short peptides, stabilizing 
folded structures such as or-helices and /5-strands which have catalytic possibilities as well as 
the association of these folded peptide subunits into larger structures, such as transmembrane 
channels. 

In analogy with modern proteins, the stability and catalytic properties of peptides will 
depend upon their ability to adopt well-defined structures. Unfortunately, short peptides 
are usually disordered in aqueous solutions and therefore, do not appear to be suitable 
for carrying out these functions. However, many of these peptides can acquire a broad 
range of well-defined secondary structures, such as a-helices, /5-strands or /5-turns, at water- 
membrane, water-oil or water-air interfaces, the precise secondary structure being determined 
by the sequence of amino acids. A crucial, common characteristic of these interfaces is that 
a nonpolar phase is adjacent to an aqueous phase. This research project has investigated 
the properties of these interfaces and its consequences on the equilibrium distribution and 
conformational and orientational ordering of small solutes and peptides at these interfaces. 

In this report, we discuss the interaction of small solutes with membrane interfaces, and 
show how the structural and electrical properties of the membrane affect the adsorption and 
conformational equilibria of the solutes. These ideas are then applied through simulations 
of small peptides at a water-hexane interface. The close proximity of the aqueous and non- 
polar phases induce the peptide to adopt conformations that are amphipathic such that the 
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hydrophilic residues are solvated by the aqueous phase and the hydrophobic residues are 
solvated by the nonpolar phase. These same principles are shown to govern the folding of 
a small, hydrophobic peptide at a water-hexane interface and the insertion of this peptide 
into the hexane phase. We then discuss the unassisted transport of ions across a simple 
glycolipid membrane. Simple solvation models, based on the different dielectric constants of 
the membrane interior and the aqueous phase, predict free energy barriers that are several 
orders of magnitude larger than inferred experimentally. Simulation results show how defor- 
mation of the membrane and transport of several solvated water molecules with the ion can 
lower this free energy barrier to a value consistent with experimental results. Finally, we 
present some preliminary results on ion channels, which are transmembrane peptide assem- 
blies which allow transport of ions across membranes with a much lower free energy barrier 
than unassisted transport. For more detailed overviews of this work, please s^e (Pohorille 
et al., 1996) and (Pohorille et al., 1999). 

Methods 

The MD simulations were carried out in the (N,V,E) or (N,p,T) ensembles. In an (N,V,E) 
ensemble, the number of atoms or molecules, N, in the system are constant and the volume, 
V, is fixed. The equations of motion for the particles in the system (Newton’s equations) are 
conservative and the total energy, E, is a constant of the motion. In an (N,p,T) ensemble, 
the volume is allowed to fluctuate around a constant fixed pressure, p, and the energy is 
allowed to fluctuate under the condition that the temperature of the system is a constant. 
For example, most simulations are carried out at temperatures close to 300 K and pressures 
close to 1 atm, which correspond to typical conditions on the contemporary earth. It is 
under these conditions that most laboratory experiments are carried out. For most systems, 
the ensemble is chosen as a matter of convenience. 

The system consisted of either one or two lamellae of water in contact with either a 
lamella of hexane or a bilayer membrane consisting of glycerol- 1-monooleate (GMO) or 1- 
palmitoyl 2-oleoyl sn-glycero 3-phosphatidylcholine (POPC), and one solute molecule. The 
actual system sizes varied, depending upon the the size of the solute molecule. 

Simulations in the (N,V,E) ensemble were: a system containing one small solute molecule, 
such as a methane derivative or small alcohol or alkane, and a lamella 500 water molecules in 
contact with a lamella of 84 hexane molecules in a simulation box with a cross sectional area 
of 24 A x 24 A; a system containing larger solutes, such as an undecamer of poly-L-leucine, 
interacting with 1380 water molecules and 409 hexane molecules in a simulation box of 42 A 
x 42 A; a water-GMO membrane system contained two lamella of 1152 water molecules, 
each, in contact with a lipid bilayer containing 72 GMO molecules in a box with a cross 
sectional area of 37 A x 37 A. Calculations carried out in the (N,p,T) ensemble used 5743 
TIP3P (Jorgensen et al., 1983) water and 64 POPC molecules in a simulation box with an 
average cross section of 51 A x 40 A. 

Periodic boundary conditions were applied in all three directions. The equation of mo- 
tions were solved using the Verlet algorithm with time steps of 2 or 2.5 fs. Calculations were 
performed at either 300 K or 310 K. In all calculations performed in the (N,V,E) systems, 
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the water was represented using the TIP4P model and the hexane molecules using the the 
OPLS (Jorgensen et al., 1984; Jorgensen and Tirado-Rives, 1988) model. The hydrocar- 
bon tail of GMO was also represented using the OPLS (Jorgensen et al., 1984; Jorgensen 
and Tirado-Rives, 1988) model and the head group potentials were derived from quantum 
mechanical calculations on glycerol. (Wilson and Pohorille, 1994) The potentials for small 
solutes were developed from quantum mechanical calculations. (Pohorille and Wilson, 1996) 
The nonbonded and the intramolecular parameters for the peptides were taken from the 
AMBER force field of Cornell. (Cornell et al., 1995) In the (N,V,E) simulations, the inter- 
molecular interactions involving water molecules and/or small, electrically neutral groups 
of the solute and hexane molecules were smoothly truncated between 7.5 and 8.0 A. For 
the (N,p,T) simulations of a water-POPC bilayer system, the head groups of the phospho- 
lipid are zwitterionic, and long-ranged forces were calculated using the Particle Mesh Ewald 
method. (Essmann et al., 1995) 

Results 

Properties of Aqueous Interfaces and Adsorption of Small Solutes: Understanding 

the structure of aqueous interfaces is the first step towards proposing physically motivated, 
implicit models of these systems. Probably the most direct characteristics of the interface 
are the mass and atomic density profiles in the direction perpendicular to the interface. The 
mass density profiles are shown in Fig. 1 and the atomic density profiles for water-octanol 
and water-POPC systems are shown in Fig. 2. The atomic density profiles for water-hexane 
and water-GMO interfaces have been published previously. (Wilson and Pohorille, 1994; 
Pohorille, 1998) In all cases, the origin was set at the equimolar surface of water, defined as 
the surface where the excess density on the side in contact with the organic phase is balanced 
by the depletion of density on the bulk water side. As can be seen from Fig. 1, the mass 
density across the interfaces is not only non-uniform, but also non-monotonic. Furthermore, 
the profiles are different for different interfaces. As will be argued below, these differences 
have direct bearing on the excess chemical potential of solutes in the interfacial region. 

Although atomic densities of the two phases in contact overlap in all four cases, it does 
not always imply interpenetration at the molecular level. In fact, simulation results have 
shown that most liquid-liquid interfaces are locally sharp, but they are broadened by spatial 
and temporal averaging over fluctuations caused by capillary waves. (Pohorille and Wilson, 
1993; Benjamin, 1997) This picture has been borne out by the results of recent scattering 
experiments on liquid-liquid interfaces. (Mitrinovic et al., 1999; Zhang et al., 1999) In the 
water-hexane system, dewetting of the nonpolar oily phase produces a small minimum in 
the mass density profile at the interface. In the water-octanol system, the interface is also 
locally sharp, but the minimum disappears due to favorable interaction between hydroxyl 
groups of water and octanol. The additional structure of the profiles for octanol is a result of 
the interfacially induced, “tail-to-tail” ordering. of octanol which persists to some extent for 
a few molecular layers. In particular, the minimum at approximately 10 A from the interface 
corresponds to the region in which tails of oppositely directed layers of octanol molecules 
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Figure 1: Mass density profiles at the water-hexane, water-octanol, water-GMO and water- 
POPC interface. The origin is defined by the equimolar surface of water, and the water lies 
to the left. 

are in contact. 

In contrast to liquid-liquid interfaces, overlaps in atomic density profiles observed for 
water-membrane system reflect penetration of water into the lipid head groups. This, in 
turn, creates a high density interfacial region, especially near very polar POPC head groups. 
Note, that density profiles for water-POPC system converge on the aqueous side to a slightly 
different value than the other three profiles because a different model of water was employed 
in the simulations. 

The non-uniform orientational distribution of molecules in the interfacial region generates 
the electric field, E z (z ), in the 2 -direction. E z (z) can be conveniently calculated as: (Wilson 
et ai, 1987) 


B,{z) = MsMzSi M (1) 

where q + (z) and <?_(z) are total charges on both sides of the plane located at height 2 and 
A is the xy cross-sectional area of the simulation box. 

The average electric fields across the interfaces are shown in Fig. 3. As expected, they 
disappear in the bulk phases but they are quite large in the interfacial regions. As was the 
case with the density profiles, the electric fields are non-monotonic and depend upon the 
nature of the interface. For the water-hexane system, E z is positive on the water side of 
the interface, but negative of the hexane side. This reflects the orientational preferences of 
water molecules in the interfacial region. The most common orientation of water molecules 
directly in contact with the organic phase is such that one 0-H bond points toward the 
organic liquid whereas the other 0-H bond is directed into the aqueous phase. 
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Figure 2: Density profiles of (a) water and octanol 0 atoms and CH n groups at the water- 
octanol interface and (b) water, POPC head group atoms and POPC tail group atoms at 
the water-POPC interface. The geometry is defined in Fig 1 


For the other three interfaces, E z is a sum of contributions from water and the head 
groups interacting with the aqueous phase. As a result, the positive peak on the aqueous 
side of water-octanol and water-GMO systems disappears, whereas the negative peak not 
only remains in place but becomes deeper. The profile of E z for water-POPC system is even 
more complicated. The minimum is shifted towards the aqueous side with respect to the 
equimolar surface of water due to the increased penetration of water into the POPC head 
groups. The additional peak on the membrane side reflects the molecular structure and 
relatively rigid orientation of the head groups. 

In general, it may be anticipated that the large interfacial electric fields will influence 
properties of solutes at interfaces. In particular, it may be expected that they will induce 
orientational preferences of small dipolar species in the interfacial region. We will show that 
this is indeed the case. 

Calculating the chemical potentials of small solutes. We are concerned with the 
calculation of the excess chemical potential, A // exc , of a solute in an interfacial environment. 
This quantity represents the quasistatic work needed to bring a solute molecule from the gas 
phase to solution and can be expressed as: 

A^X e xc fA f^idealj (2) 

where fi and fj.\deai are the chemical potentials of the solute in solution and in ideal gas, 
respectively. In interfacial systems, A// exc is a function of the ^-coordinate, perpendicular to 
the interface. However, since this dependence has no influence on the derivations presented 
below, we will omit it to simplify notation. 

A very efficient method for estimating A/j, exc of small solutes is provided by the potential 
distribution theorem. (Widom, 1963; Widom, 1982) In this method, only the system without 


7 




Wilson, Michael A. 
NCC2-772 



Figure 3: ^-component of the average electric field at aqueous interfaces. The geometry is 
defined in Fig 1 


the solute needs to be simulated. Then, advantage is taken of the relationship between 
the excess chemical potential of the solute inserted into the system and the corresponding 
solute-solvent potential energy, U so i , 

A// e rc — kgT In < exp( Usol/kgT ') > 0 ) (3) 

where kg is the Boltzmann constant, T is temperature and < ... > 0 signifies a thermal 
average over insertions into configurations of the system unperturbed by the solute. Here, 
we assume that the solute is rigid, a good approximation for the molecules considered in 
this paper. Furthermore, Eqn. 3 is exactly correct only for the NVT ensemble. In the NVE 
ensemble, a correction arising from temperature fluctuations should be included. (Frenkel 
and Smit, 1986) This correction becomes important near the critical temperature, but under 
conditions of the present simulations is expected to be small and, therefore, has been omitted. 

Eqn. 3 simplifies greatly if the solute is a hard sphere (a cavity). Then, exp(-U sol /kgT) 
is either 1 if the inserted hard sphere does not overlap with van der Waals spheres around 
solvent atoms or 0 otherwise and 

A/ierc = -k B T\n[P cav (R)], (4) 

where P cav {R) is the probability of inserting a cavity of radius greater than or equal to R 
into the system. One efficient method of calculating P ca v{R ) has been proposed by Pohorille 
and Wilson, (Pohorille and Wilson, 1996) and the same method has been used in this work. 

Once P C av(R) is known it can also be used as the probability density for inserting solute 
molecules. Then, Eqn. 3 takes the form: 
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A/^eXC — 



Pcav(R) < exp(-U so i/k B T ) > R dR . 


(5) 


where < ... > B represents a statistical overage over insertions of the solute into cavities with 
radii between R and R + dR. 

Since Pcav(R) is a rapidly decreasing function of R the efficiency of such a particle inser- 
tion method would be poor. Most cavity sites sampled from this distribution are too small to 
accommodate the solute, yielding very large U so i, and, consequently, contribute practically 
nothing to the statistical average in Eqn. 5. However, the efficiency can be greatly improved 
by introducing a probability distribution, P(,, biased by a weighting function, w(R), which 
increases with R: 


V 

Pb(R) = Pcav(RMR). 


Then, Ap exc can be expressed as: 


(6) 


Ap exc 


-k B T In 


'j£ 


Rmin 


ncg) 

w(R) 


< exp(-U so i/k B T ) > R dR \ , 


(7) 


where is the smallest cavity size considered a suitable insertion site. 

To calculate Ap exc as a function of 2 the system is divided into N layers of width A 2 and 
the averaging in Eqn. 6 is carried out in each layer separately. 

Calculations of p exc (z) start from configurations generated in molecular dynamics simu- 
lations of the neat interfacial system (without the solute). For each of these configurations a 
large number of solute insertions are made, governed by the probability distribution Pb(R), 
and U so i is evaluated at each insertion site. These potential energies are further used to 
estimate the statistical average in Eqn. 7. 

A major advantage of the particle insertion method is that it may yield p exc of several 
different solutes at many different locations along 2 from a single simulation of the system 
in the absence of these solutes. This makes the method both highly efficient and accurate, a 
rare combination in computer simulations. The method, however, also has a major limitation 
— it is applicable only to atomic and small molecular solutes that can fit into voids that 
transiently form in the system. If the solute is too large, the number of sterically allowed 
insertions would be small and the accuracy of estimating p exc would suffer. 

Excess chemical potentials of small solutes: In implicit models, the excess chemical 
potential of a solute molecule is commonly represented as a sum of two contributions — 
electrostatic and non-electrostatic — which are evaluated separately. The non-electrostatic 
term is usually obtained from empirical formulas for the reversible work required to create 
a cavity of appropriate size in the solvent, and for solute-solvent van der Waals interactions. 
In explicit models, the total excess chemical potential is calculated at once. Nevertheless, it 
is still possible to separate electrostatic and non-electrostatic terms by performing properly 
designed simulations. This is exactly the approach that we will follow. This approach not 
only may guide future developments of interfacial implicit models but also offers valuable 
insight into the origin of the observed profiles of Ap exc . 
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Figure 4 : Excess chemical potentials of methane and its fluorinated derivatives across the 
(a) water-hexane, (b) water-octanol, (c) water-GMO, and (d) water-POPC interfaces. The 
geometry has been defined in Fig. 1. 


First, we discuss the profiles of total A (i exc , shown in Fig. 4. For the three dipolar 
fluoromethanes (CH3F, CH2F2 and CHF3) the profiles exhibit interfacial minima, indicating 
that these molecules tend to accumulate at the interface. The profiles, however, are not 
identical for all four interfacial systems. At the water-hexane interface, the minima are deeper 
(by approximately 1 kcal/mol) but not as broad as those at the water-GMO interface. The 
profiles in the water-octanol system reveal a lot of structure which is related to the extended 
layering of octanol in the interfacial region. However, due to difficulties in fully equilibrating 
octanol, it is uncertain whether the small features in these profiles away from the interface are 
statistically significant. For the water-POPC system, A fi exc exhibits not only a minimum, 
which is located on the hydrocarbon side of the head group region, but also a maximum 
on the water side of the interface. This feature is not present in the other three interfacial 
systems. 

The profiles for nonpolar molecules, CH 4 and CF 4 , are qualitatively different from the 
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Figure 5: Excess chemical potential of cavity formation across the water-hexane, water- 
octanol, water- G MO, and water-POPC interfaces. The geometry has been defined in Fig. 1. 

polar species. They exhibit a weak minimum at the water-hexane interface and no interfacial 
minima in the other three systems. A/i ezc decreases monotonically or nearly monotonically 
across the water-octanol and water-GMO interfaces and exhibits a maximum at the water- 
POPC interface. A weak minimum in octanol separated from the interface by approximately 
one molecular layer reflects, again, interfacially induced ordering of octanol. In summary, in 
contrast to their polar analogs, nonpolar solutes concentrate in the organic phase or in the 
middle of the membrane bilayer. 

The reversible work of cavity formation, fj, c av , which provides the main contribution to 
the non-electrostatic term, can be calculated from Eqn. 4 and the full non-electrostatic term 
can be obtained from Eqn. 3 by inserting into the solvent a hypothetical, nonpolar solute 
which has the same van der Waals parameters as the solute of interest but no distribution of 
partial charges. Then, the electrostatic term can be obtained by subtracting A fi exc of this 
hypothetical solute from A fi exc of the solute with the correct charge distribution. 

In Fig. 5, we show fi cav (z) for cavities of radius 2 A across the four interfaces studied in 
this work. The most salient feature of these profiles is that they are non-monotonic in the 
interfacial region and change differently for different interfaces. Comparison between Figs. 1 
and 5 reveals similarities between fi C av and the mass density profiles for the corresponding 
interfaces. These similarities are not surprising. It is expected that it should be easier to 
insert a hard sphere in low mass density interfacial regions than in high mass density regions 
because they contain more free volume. The water-hexane interface serves as an example of 
such a region. In the water-octanol system, the minimum in fi ca v(z) is not located exactly at 
the interface, but instead in the low density tail region between the first and second layer of 
octanol. In contrast, when mass density reaches a maximum at the interface, as is the case 
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for the water-POPC system, we expect to find an interfacial maximum in fi ca v(z)- It should 
be, however, kept in mind that this is only a qualitative argument and does not capture 
certain important structural features of solvents. For example, it has been shown that it is 
more difficult to insert a cavity of atomic or small molecular size into water than into other 
liquids even though they may have less free volume. (Pohorille and Pratt, 1990; Pratt and 
Pohorille, 1992). This is a reflection of the hydrophobic effect and it can be, observed in 
Fig. 5 by comparing fi cav on both sides of the interfaces and noting that it is always larger 
on the water side. 

The inclusion of van der Waals interactions in A/j, exc shifts the profiles of excess chemical 
potential towards lower values, but has little influence on their overall shape. Thus, in 
discussing the behavior of small, nonpolar species at different aqueous interfaces, we can be 
guided by p cav {z). v 

For other nonpolar molecules, profiles of Ap exc can take somewhat different shapes. One 
possible reason for such differences is that the reversible work of cavity formation may scale 
differently with cavity size in different regions along the z-direction. (Marrink et al., 1996) 
Additionally, the ordering of phospholipid chains influences the shape of cavities. This is 
particularly evident in the region adjacent to the head groups. In this region finding a non- 
spherical cavity aligned with the lipid chains is considerably more probable than finding a 
spherical cavity with the same volume. (Bassolino-Klimas et al., 1993; Marrink et al., 1996) 

The profiles for the electrostatic contributions, A/i e ;, can be obtained by subtracting 
Afi exc for hypothetical solutes without partial atomic charges from A fj, exc shown in Fig. 4. 
This approach was taken in the previous analysis of solutes at the water-hexane and water- 
GMO interfaces. (Pohorille and Wilson, 1996) This analysis, however, can be further clarified 
if, instead of the molecular solute, we consider a model solute consisting of a point dipole 
in the center of a cavity. This model is sufficiently simple that several useful results can 
be derived analytically. Simultaneously, as has been already shown (Pohorille and Wilson, 
1996), A /u exc (z) for fluoromethanes can be faithfully reproduced by such a model if the cavity 
radius is properly chosen. 

For the dipolar model, the solute-solvent electrostatic energy, U e i, is given as: 

U e i = — p • £ = —p£cos0, (8) 

where p is the dipole moment of the solute and £ is the electric field created by the solvent and 
acting on the dipole and 9 is the angle between these two vectors. Then, Afi e i, statistically 
averaged over dipolar orientations, is expressed as: 


A W = -ksTlu < >„= -k B T In < >0 . (9) 

Here < > 0 represents a statistical average over configurations of the reference system 

(without the dipole). 

In Fig. 6 we show interfacial profiles of Ap t i for a 2.0 D dipole in a cavity the radius of 
which is 2 A. The most striking feature of these profiles is that, within accuracy of the data, 
they are all quite similar. A// e ; is constant in water, then increases approximately linearly in 
the interfacial region over a distance of 10 A, and finally reaches a plateau of approximately 
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Figure 6: Electrostatic contribution to the excess chemical potential of a dipole of 2.0 D in a 
cavity of radius 2.0 A across the water-hexane, water-octanol, water-GMO, and water-POPC 
interfaces. The geometry has been defined in Fig. 1. 


zero in the bulk organic phase or inside the membrane bilayer. The shape or the width of 
the interfacially changing portion of the profile does not appear to depend significantly on 
the structure, polarity or ordering of chemical groups directly in contact with water. 

Another issue of particular relevance in the context of developing successful implicit mod- 
els is whether the linear response approximation is accurate at the interface. This requires 
that the second order perturbation theory for Ap e / is satisfactory. In this approximation, 
A n e i is expressed as: 


k „T 

A fid =< U e , > 0 — —[< Ul > 0 - < U el >2] (10) 

Since < U e i >o averaged over dipolar orientations vanishes this expression simplifies to: 

A^ el (z) = < S(z) 2 > 0 . (11) 

In the last formula we explicitly included the dependence on z to stress that A p e / changes 
with the distance from the interface. 

Eqn. 1 1 implies that Afi e i should be proportional to p 2 and, therefore, can be considered as 
a molecular-level equivalent of the well known Onsager formula obtained from the continuum 
model of bulk solvents. (Onsager, 1936). The dependence of A p e / on p 2 in bulk aqueous 
solution and in the interfacial regions is shown in Fig. 7. As can be seen, second order 
perturbation theory is quite accurate in water, but is less satisfactory at the water-hexane 
and water-octanol interfaces. Surprisingly, the linear relation between A p e / and p 2 appears 
to be well preserved at the water-POPC interface, despite large electric fields acting in this 
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Figure 7: Electrostatic chemical potentials of dipoles in a cavity of radius 2 A as a function 
of the square of the dipole moment in bulk water and at the water-hexane, water-octanol 
and water-POPC interfaces. 


region. This may be a result of fortuitous cancelation of errors. Alternatively, it is possible 
that the accuracy of second order perturbation theory improves if long range effects are 
fully taken into account, as long-range interactions were only included in the water-POPC 
simulation. 

Finally, we consider the orientational preferences of solutes at the interface induced by 
electric fields generated by the solvent. To this end, we examine point dipoles placed in a 
cavity and pointing parallel, antiparallel or perpendicular to the excess electric field, E z , 
acting along the z-axis. The corresponding chemical potentials are A/x t , A/4 and A/xj., 
respectively. In the second order perturbation approximation they are expressed as: 

A/'t = ~P< £|| >0 “2 E \\ >° - < E \\ >0] 


A tH = P< £|| >0 E \\ >0 - < £|| >ol 


A fix = 


2k B T 


< El > 0 


( 12 ) 


where E\^ and E± are components of the instantaneous electric field acting on the dipole in 
the directions parallel and perpendicular to E z . 

If only the first order terms are considered, A (i exc of the dipole oriented parallel to E z 
should decrease proportionally to the excess electric field at location 2, A/x erc of the dipole 
oriented antiparallel should increase by the same value and A// exc of the dipole oriented 
perpendicular should be unaffected by the electric field. In other words, if perturbation 
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Figure 8: Excess chemical potentials of dipoles of 2.08 D in cavities of radius 2.0 A oriented 
parallel, antiparallel and perpendicular to the interface across the water-POPC interface. 
The average of the excess chemical potentials of the parallel and antiparallel dipoles is also 
shown. The geometry has been defined in Fig. 1. 


theory to first order is accurate, then the excess chemical potential of the perpendicular 
dipole should simply be the cavity free energy and the excess chemical potentials of the 
parallel and antiparallel dipoles should lie symmetrically below and above this value. 

The corresponding numerical results for the water-POPC system are shown in Fig. 8. 
In the most general case, the electric field can change sign in the interfacial region, as we 
observed at the water-hexane interface (see Fig. 3), and the above notation would have to 
be generalized for orientations of dipoles that are fixed with respect to the interface. Since 
the electric fields at the water-POPC are of a single sign, to within our numerical accuracy, 
we retain the simpler notation. Although the lack of orientational averaging limits statistics 
and yields somewhat noisy profiles, certain features of these profiles are clear. First, dipoles 
at the interface tend to align themselves with the electric field. Second, the first order ap- 
proximation captures the qualitative features of the profiles, and A^ exc of the dipole oriented 
parallel to the electric field lies below A/i exc of the dipole in the antiparallel orientation in 
most of the interfacial region. The average of A fi exc for the parallel and antiparallel dipoles 
has the same shape as A fi exc for the perpendicular dipole, but is consistently higher. This 
means that the second order term makes a significant contribution to the total chemical 
potential. It also indicates that fluctuations of the electric field are larger in the plane of the 
interface than in the perpendicular direction. 

Organization of Peptides at Membrane Interfaces. We had previously studied the 
interfacial behavior of very short, model peptides and longer peptides composed of both polar 
and nonpolar amino acids. (Pohorille and Wilson, 1994a; Pohorille and Wilson, 1994b) All 


15 


Wilson, Michael A. 
NCC2-772 


these peptides exhibit interfacial activity, i.e. they tend to accumulate at aqueous interfaces. 
Recently, we extended these studies to the undecamer of poly-L-leucine at the water-hexane 
interface. Since this peptide consists only of nonpolar residues, it is not expected to be 
interfacially active but, instead, to partition into the nonpolar phase. Thus, our calculations 
consider a simple model of peptide insertion into the interior of a membrane, a phenomenon 
fundamental to the transduction of signals and the transport of material between the cell 
interior and the environment. 

It is expected that most peptides composed of mixtures of polar and nonpolar residues 
will adopt amphiphilic, interfacially active structures. Whether these structures are ordered 
or not depends on the sequence of amino acids. In contrast, nonpolar peptides are expected 
to partition from water to a nonpolar medium. In aqueous solution, such peptides are 
disordered, whereas in nonpolar media, they most likely exist as a-helices. They must, v 
therefore, fold during transfer across the interface. This is clearly relevant to the formation 
and insertion of transmembrane helices. A proposed scheme for this process is that integral 
peptides first accumulate at the water-bilayer interface and form the helix. (Jacobs and 
White, 1989) Once folded, the peptides penetrate the bilayer via the nonpolar regions of the 
a-helices. Alternative schemes, in which penetration precedes, or is coupled to, folding can 
also be envisioned. 

We investigated these possibilities in the example of poly-L-leucine. First, we tested the 
predictions about its stability in water and hexane. When the undecamer was placed in water 
as a /^-strand, it rapidly unfolded into a random coil — a family of disordered conformations 
observed during 6.2 ns of a molecular dynamics trajectory. The final structure was used 
in subsequent simulations as a model for the random coil. In contrast, when the peptide 
was place in hexane in a /3-strand or a random coil conformation, it refolded into an a- 
helical structure within 3.0 and 3.7 ns, respectively. In both cases, the folding proceeded 
sequentially from the C-terminus. 

To simulate poly-L-leucine at the water-hexane interface, the peptide, in a random coil 
conformation, was placed in water such that its center of mass was ca. 11 A from the 
interface. Initially, the undecamer moved rapidly towards the interface, partitioning to 
the interface within 1 ns. This indicates that a strong average force attracts a nonpolar 
peptide to the nonpolar phase. Once the peptide reached the interface, it started slowly 
folding into an a-helix and, simultaneously, partitioning into hexane. These two processes 
appear to be closely coupled: some penetration into hexane is needed to initiate folding. 
Then, as the helical structure started forming, some hydrophilic atoms of the backbone 
became shielded from the aqueous environment. This made the peptide structure more 
hydrophobic, allowing for greater penetration into the nonpolar phase and promoting further 
folding. The correlation between the evolution of the folding of the backbone and the center- 
of-mass (COM) of the peptide are show in Figure 9. In the present study, folding was not 
sequential. First, six consecutive residues, from the N-terminus, formed a helical fragment. 
Subsequently, four residues, from the C-terminus, folded into a helix. Here, the term “helix” 
denotes the portion of the Ramachandran map, for which —180 < <f> < 0 and —90 < ip < 0. 
After 34 ns, the completely formed helix was immersed in hexane, just next to the interface. 

Folding of the peptide proceeded through an intermediate, a 3i 0 -helix. During the sim- 
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Figure 9: Position of the peptide center-of-mass (dotted line) and the deviation of the peptide 
backbone from an a-helix (solid) line as a function of time. 


ulations, this helix persisted for only 2 ns before converting to the a-helix. The two helices 
differ in the pattern of hydrogen bonding along the peptide backbone. In an a-helix, each 
hydrogen bond involves residues separated by three other residues, whereas in a 3io-helix 
the participating residues are separated by only two other residues. Occasional formation of 
the 3io-helix was also observed for the folded peptide. This indicates that this helix not only 
mediates folding but remains in equilibrium with the a-helix once this process is completed. 
A very similar conclusions was reached from experimental studies of alanine-based peptides 
in aqueous solution. 

The orientation of the peptide, defined by the direction of the end-to-end vector, is nearly 
parallel to that interface. However, orientations nearly perpendicular to the interface with 
the N-terminus buried in hexane are also probable. In contrast, perpendicular orientations 
with the C-terminus in hexane are highly unfavorable since they require dehydration of 
the three carbonyl group at this end of the peptide which are not involved in intramolecular 
hydrogen bonding. Most transmembrane proteins exhibit similar preferences and incorporate 
into the membrane through the N-terminus. 

Since nonpolar peptides appear to adopt helical structures in nonpolar environments, in- 
dependently of their specific sequence, we propose that inserting such peptides into nonpolar 
liquids and membranes requires simultaneous folding and penetration. Although results from 
a single trajectory are insufficient to draw general conclusions about the mechanism of fold- 
ing, we note that nonsequential pathways are, at least, probable. Determining whether they 
are also the most preferred, and how equilibrium free energy and kinetic effects contribute 
to the selection of folding pathways, requires further studies. 

Transport of ions across membrane interfaces. Simple ions, such as Na + , K + , Cl - 
and Ca 2+ , play an essential role in a wide variety of cellular processes, such as bioener- 
getics, intercellular signaling and the catalysis of chemical reactions. To enter cells, ions 
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must permeate their lipid membranes. There is a large free energy of activation associated 
with transferring a charged species from a polar environment to the nonpolar interior of a 
membrane. Thus, a membrane forms a permeability barrier to ions, making them excel- 
lent structures for mediating many cellular regulatory processes. In contemporary cells, ion 
transport is precisely controlled with the assistance of transmembrane ion channels, carriers 
and pumps located in membranes. However, to fully understand the transport of charges 
across membranes, it is necessary to begin with passive ion transport. Such uncontrolled 
transport is usually unproductive, and sometimes harmful to cells, but useful in artificial, 
membrane-bounded structures such as liposomes. 

The permeation rate of ions across membranes can be estimated using a continuum 
dielectric model of a water-membrane system. In this model, both water and membrane 
are represented as homogeneous, isotropic media, characterized by dielectric constants e w 
and Cb, respectively, and separated by a sharp planar boundary. If the ion is represented 
as a point charge q located at the center of a cavity of radius a, the change in the excess 
chemical potential associated with the transfer of the ion from bulk water to the center of 
the membrane (the free energy barrier), A//J x s c , is expressed in this model as (Parsegian, 
1969; Neumcke and Lauger, 1969): 




C6 




^logPM 


(13) 


where w m is the width of the membrane. For an ion with the effective radius of Na + , A pj x s c is 
approximately equal to 47 kcal/mol. This yields a membrane permeability of 10“ 29 cm s -1 , 
approximately 17 orders of magnitude lower than the measured permeability of Na + (Wilson 
and Pohorille, 1996). This discrepancy is so large that it cannot be easily explained by 
uncertainties in the parameters of equation (13). Thus, the continuum dielectric model 
fails to accurately represent the energetics of ion transfer across biological membranes. The 
question is, why does it fail? 

To help answer this question, a detailed, molecular simulation of the unassisted ion trans- 
port of Na + and Cl across a GMO bilayer was performed using molecular dynamics (Wilson 
and Pohorille, 1996). Due to the large free energy barrier, however, the transfer does not 
occur on the timescales accessible to computer simulations. Therefore, a series of molecular 
dynamics trajectories were generated in which the ion was restricted to overlapping regions 
along the interface normal (^-direction). 

The simulations revealed a picture of ion permeation that is in sharp contrast with the 
continuum dielectric model. As the ion moves across the water-membrane interface into 
the bilayer, the membrane surface does not remain approximately planar. Instead, a local 
deformation is formed in which water molecules and polar head groups (normally restricted 
to the surface of the membrane) follow the ion into the nonpolar interior of the bilayer. Once 
the ion crosses the midplane of the membrane, the deformation on the incoming side relaxes 
and simultaneously, a similar deformation forms on the outgoing side. Thus, during the entire 
transfer process, the ion remains partially solvated by both the polar head groups and water 
molecules. The key feature of this molecular description of the ion transfer process is that 
the ion is never fully solvated by the nonpolar hydrocarbon tails. Thus, the calculated A/zJ x s c 
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is markedly lower than the barrier predicted from the continuum model. For Na + , A/zJ x ^. 
was estimated at 27 kcal/mol. The full profile of A fi exc across the membrane is shown in 
Figure 10. The barrier for Cl - was approximately 2 kcal/mol lower. To obtain this estimate 
it was necessary to correct for the absence of permanent and induced partial charges on the 
atoms of the GMO tails. The lack of charges would yield ej of the membrane interior equal 
to 1 instead of 2, as expected for real membranes (Ohki, 1968; Dilger et al. , 19§2). 



Figure 10: Free energies of transferring Na+ across the water GMO interface. The center of 
the membrane is at z = 0. 

Once A/i^ x s c is known, the permeability coefficient can be estimated by expanding 
A// e xc(^) in the high barrier limit (Wilson and Pohorille, 1996): 

Pmembr = Db exp ^ (14) 

where a = (1/2) | <9 2 A/i eX c {z)jdz 2 | and Dt, is the diffusion coefficient in the z-direction at 
the midplane of the bilayer. 

The calculated P mem br for Na + is approximately 10 -15 cm s -1 . This value is slightly lower 
than most experimental estimates for phospholipid bilayers (Hauser et al., 1973; Nozaki and 
Tanford, 1981; Deamer and Volkov, 1995). Considering, however, that the error in the 
calculated A/xJ x s c is approximately 4 kcal/mol, corresponding to an uncertainty in P membr of 
3 orders of magnitude, the agreement between the simulated and measured permeabilities is 
satisfactory. The permeability coefficient for Cl - , computed from the molecular dynamics 
simulations, is approximately 2 orders of magnitude higher than that for Na+ , again in 
good agreement with experiment (Hauser et al, 1973; Paula et al, 1998). The difference in 
permeability between Na + and Cl - most likely reflects the difference in fj. exc of these two ions 
in water. This, in turn, is associated with the larger ionic radius of Cl - , properly adjusted 
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for a different orientation of water molecules around positive and negative ions (Rashin and 
Honig, 1985). As a result, partial desolvation accompanying the transfer of an ion from 
water to the membrane requires less free energy for Cl~ than for Na + . 

These simulations were performed only for one type of the membrane. Therefore, there 
is no direct data on the effect of different head groups and membrane widths on ionic per- 
meability. It can be expected, however, that for sufficiently wide membranes, the formation 
of deep deformations, extending to the midplane of the bilayer, would become highly unfa- 
vorable. Then, the ion would have to undergo almost complete desolvation near the center 
of the bilayer. In that region, the change in A^J X S C may become similar to that predicted 
from the continuum dielectric model. 

Conclusions v 

Computer simulation of several aqueous interfaces, from the simple water-hexane interface 
to the more densely packed water-POPC interface have show that the adsorption of small 
solutes at these interfaces can be understood as a balance between electrostatic effects and 
non-electrostatic terms. For polar solutes, the free energy arising from electrostatic effects de- 
creases from the nonpolar phase to the water. On the other hand, the non-electrostatic term 
increases from the nonpolar phase to water. We find that the electrostatic term decreases 
sufficiently strongly to give rise to a free energy minimum at water-oil and water-membrane 
interfaces, indicating that small, polar solutes are adsorbed at these interfaces. 

Separating the free energy into its electrostatic and non-electrostatic parts also highlights 
two important points. First, the variation in the electrostatic term in the interfacial region 
is similar for all the interfaces we have investigated. This suggests that it might be possible 
to build simple implicit theories to account for this term. Second, the non-electrostatic term 
is a sensitive function of the solvent density in the interfacial region. This term shows a free 
energy minimum at the water-hexane interface, where the hydrophobic phase is directly in 
contact with the aqueous phase. But this minimum disappears for the more densely packed 
amphiphilic, bilayer systems. This is a significant problem because it is not clear that it will 
be possible to extend theories of cavity formation, which have been well characterized for 
bulk phases, into the interfacial region. 

The electric fields and free energy required to form cavities at the interface clearly influ- 
ence the adsorption and and orientational and conformational equilibria of small solutes in 
the interfacial region. For larger molecules, such as peptides, the situation is more complex. 
If peptides consist of nonpolar residues only, they become inserted into the nonpolar phase. 
As demonstrated by the example of the leucine undecamer, such peptides fold into an c*-helix 
as they partition into the nonpolar medium. The folding proceeds through an intermediate, 
called the 3io-helix, which remains in equilibrium with the a-helix. This process represents 
a simple, protobiologically relevant example of environmentally-mediated self-organization 
of biological molecules. Once in the nonpolar environment, the peptides can readily change 
their orientation with respect to the interface from parallel to perpendicular, especially in re- 
sponse to local electric fields. The ability of nonpolar peptides to modify both the structure 
and orientation with changing external conditions may have provided a simple mechanism 


20 


Wilson, Michael A. 
NCC2-772 


of transmitting signals from the environment to the interior of a protocell. 

While the free energy of cavity formation and the interfacial electric fields determine the 
permeability small polar solutes, they have a much smaller effect on the unassisted transport 
of ions across a membrane. In this case, it is the flexibility of the bilayer that determines 
the mechanism of ion transport across asymmetric thinning defects in the membrane. In 
this mechanism, the head groups on the incoming side of the bilayer follow the movement 
of the ion into the membrane interior by tilting inwards. The resulting defect is filled with 
water. Once the ion crosses the mid-plane of the bilayer, the defect on the incoming side 
of the membrane disappears and a similar deformation is formed on the outgoing side. 
During the transfer, the ion remains solvated by both head groups and water molecules. 
Even though the calculated ionic permeabilities are subject to considerable uncertainties, 
this mechanism qualitatively accounts for a large increase in these permeabilities (by 14-17 
orders of magnitude), compared to the estimates from the rigid membrane model. Our results 
imply that the permeability of thin membranes to ions should change with the bilayer width 
more rapidly than expected from the continuum dielectric model, also in agreement with 
experimental results (Deamer and Volkov, 1995). This offers, for example, the possibility 
of regulating release of charged therapeutic agents from liposomes. The calculations also 
illustrate the broader point hat polar groups and/or molecules in the membrane can readily 
reorganize to stabilize a charge inside the bilayer. It is likely that this effect should be 
considered, for example, in the studies of the mechanism by which peptides with charged 
residues are incorporated into lipid bilayers. 

Future Directions 

One motif that is both appealing as a method of obtaining new functions from existing 
peptides and exists in many actual protein assemblies such as the superfamily of ligand- 
gated ion channels, is the aggregation of several peptides together to form a transmembrane 
assembly. Such an assembly could have functioned as a simple signaling system or as an ion 
transporter. We are investigating the free energetics and kinetics of the dimerization and 
insertion of the designed peptide, (LSLLLSL) 3 , which forms tetrameric ion channels in the 
presence of an applied external voltage. Owing to the importance of proton transport across 
membranes throughout all living, cellular systems, we have initiated studies of a simple, 
proton-conducting, transmembrane channel to understand possible mechanisms of transport 
and gating. 
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